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Excited state soliton stars are studied numerically for the first time. The stability of spherically 
symmetric S-branch excited state oscillatons under radial perturbations is investigated using a f D 
code. We find that these stars are inherently unstable either migrating to the ground state or 
collapsing to black holes. Higher excited state configurations are observed to cascade through inter- 
mediate excited states during their migration to the ground state. This is similar to excited state 
boson stars |f|]. Ground state oscillatons are then studied in full 3D numerical relativity. Finding 
the appropriate gauge condition for the dynamic oscillatons is much more challenging than in the 
case of boson stars. Different slicing conditions are explored, and a customized gauge condition 
that approximates polar slicing in spherical symmetry is implemented. Comparisons with fD re- 
sults and convergence tests are performed. The behavior of these stars under small axisymmetric 
perturbations is studied and gravitational waveforms are extracted. We find that the gravitational 
waves damp out on a short timescale, enabling us to obtain the complete waveform. This work is a 
starting point for the evolution of real scalar field systems with arbitrary symmetries. 
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I. INTRODUCTION 

Real scalar fields play an important role in many mod- 
els in particle physics and cosmology. The axion is de- 
scribed by a real field as is the inflaton. Scalar particles 
are important dark matter candidates. Models of real 
scalar field dark matter have given good fits for the veloc- 
ity profile of spiral galaxies [4]. Recently there has been 
widely publicized observational evidence by NASA for 
the existence of dark matter with both the Hubble tele- 
scope and the Chandra X-ray observatory ^4|. These 
particles could come together through some kind of Jeans 
instability mechanism to form gravitationally bounded 
objects such as oscillating soliton stars (also called oscil- 
latons) (Bl. 

Soliton stars possess a spectrum of states: ground state 
and excited states. The ground state scalar field has no 
nodes, a first excited state has one node and so on. The 
stability of spherically symmetric ground state oscillatons 
has been studied numerically [1, 0I- In this paper we 
extend this study to the case of excited states. Studying 
the stability of excited states is important because they 
may be intermediate states during the formation of these 
stars. 

In this work we also begin an exploratory investigation 
of ground state soliton stars in three-dimensional numer- 
ical relativity with the goal of performing stable evolu- 
tions and extracting gravitational waveforms. Numeri- 
cal relativity has made major progress in the past few 
years. It is now possible to evolve binary black hole in- 



spirals stably for many orbits [8|, |9| . Gravitational wave- 
forms from black hole simulations have been extracted 
and match accurately against post-Newtonian estimates 
(l0| . Gravitational wave detectors have reached their de- 
sign sensitivity and waveform templates are required to 
analyze data produced by the science runs Gen- 
erating numerical waveforms from other compact objects 
predicted by general relativity such as soliton stars (com- 
prised from real scalar fields) and boson stars (made up 
from complex scalar fields) is a timely and necessary ef- 
fort. Gravitational waveforms from two colliding boson 
stars |T^ . [1^ . ITU as well as from single distorted boson 
stars [lal have been extracted. To our knowledge this is 
the first soliton star study in 3D numerical relativity. 

Numerous authors have studied soliton stars numer- 
ically in spherical symmetry 0, [13, [3 • Mathemati- 
cally, oscillating soliton stars are nonsingular solutions to 
the Einstein-Klein-Gordon (EKG) equations represented 
by a massive real scalar field for which both the metric 
and the scalar field are periodic in time Q. The ab- 
sence of static equilibrium configurations contrasts oscil- 
latons with the case of boson stars, which are complex 
scalar field solutions to the EKG equations that have 
equilibrium configurations characterized by static metric 
components (the scalar field has an ^{r,t) = (/<(r)e*"* 
time dependence but the metric and energy itself are 
time independent), and makes their evolution more chal- 
lenging numerically. Similar to boson stars and neu- 
tron stars, the mass profile of ground state soliton stars 
has a stable (^-branch) and an unstable branch {U- 
branch) Q. The mass profile has a maximum at 
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A^cgroundst. = 0.607M|j/m Corresponding to a central 
density 0i(O) = 0.48 (here m is the mass of the scalar 
field.) The configurations to the left (lower central den- 
sity) of this critical mass are stable, while those to the 
right (greater central density) are unstable. 

In this paper we study the stability of excited state os- 
cillatons under radial perturbations and find that all ex- 
cited configurations are inherently unstable. They either 
migrate to the ground state or collapse to black holes. 
We present simulations for both scenarios for configura- 
tions in the first excited state. We find the mass profile 
of first excited state configurations as a function of cen- 
tral density (f>i{0) (see Fig. [TJ) It has a maximum at 
a similiar critical density as for ground state stars, but 
in this case it corresponds to a higher critical mass of 
Mc 1st excited St. = 1.33M|[/to. We evolve S'-branch one- 
node stars when no explicit perturbation is applied other 
than that induced by the numerical grid (we covered the 
95% mass radius of the star with the same number of grid 
points ~ 187 to make meaningful comparisons) and find 
that stars that are unable to migrate to the ground state 
collapse to black holes. The migration can occur only if 
the star can lose enough mass through scalar radiation 
to become an 5-branch ground state configuration with 
M < Mc ground St. = 0.607Afp[/m. A first excited state 
star of mass M < 0.84A/|[/m succeeds in losing its ex- 
cess mass and migrates to the ground state. We study the 
relative times of collapse to black holes of different config- 
urations for S'-branch stars with mass M > O.SAMpJm. 
The collapse times that we tabulate are approximate be- 
cause polar slicing does not penetrate black hole horizons. 
We find that this time decreases as the central density in- 
creases. We also present a simulation of the migration of 
a one-node star to the ground state and find the ground 
state configuration that it migrated to. We then simulate 
the migration of a 5-node configuration to the ground 
state. During the migration process the star cascades 
through a superposition of lower excited state configura- 
tions. We are able to determine an approximate 4-node 
intermediate state and then follow the evolution until the 
star settles into a ground state configuration. The sim- 
ulations are performed using the polar slicing condition 
Kg + = 0, which is a natural choice for spherically 
symmetric spacetimes [l9j . 

In the second part of the paper we focus on 3D simu- 
lations of ground state soliton stars. We use the Cactus 
Computational Toolkit |20|] with the scalar field evolution 
of Guzman [2l| and the BSSN implementation [l^, of 
the Einstein equations. Soliton stars retain some of the 
properties of boson star systems that made them a good 
testbed for numerical relativity in that the space-time is 
singularity free and the star has a smooth outer bound- 
ary. However, the evolutions are challenging because the 
system is highly dynamic, i.e., the metric functions, ex- 
trinsic curvature components, and gauge functions are 
always rapidly oscillating in time. In particular, find- 
ing an appropriate slicing condition and implementing it 
accurately to obtain stable evolutions is a challenge for 



such dynamic spacetimes. We introduce a customized 
slicing condition based on the features of an eigenstate 
of a stable branch oscillaton in spherical symmetry and 
we refer to it as truncated Fourier slicing. Instead of en- 
forcing that K be zero as in the case of dynamic boson 
star evolutions 16], for this system we drive K to the 
time dependent iiT^j'"'*''-' (i), which is the trace of extrinsic 
curvature appropriate for a spherically symmetric oscilla- 
ton eigenstate expanded to finite order in a Fourier series. 
While this condition is based on an eigenstate, we find 
that it is effective in simulations with small nonspheri- 
cal perturbations applied to the star. With this choice 
of gauge we are able to reproduce results from ID codes 
within the 3D context, and we are able to perform stable 
simulations of sufficient duration to extract gravitational 
waveforms. 

We then study soliton stars under small nonradial per- 
turbations and the emitted gravitational waves. We use 
two types of perturbations. The first type was previ- 
ously applied to boson stars by Guzman [2l|. It con- 
sists of an imposed asymmetry in the grid resolution 
Aa; — Ay ^ Az, while choosing the number of points 
nx — n-y Uz such that the distance from the origin to 
the boundary of the grid is kept the same. In this case the 
resolution itself is the perturbation. Consequently, the 
amplitude of the £ = 2,m = Zerilli waveform is ex- 
pected to be zero in the limit of infinite resolution and is 
observed numerically to converge away with resolution at 
approximately second order. We next study an oscillaton 
under a perturbation proportional to the i = 2,m — 
spherical harmonic that perturbs the metric of the eigen- 
state nonradially. This perturbation is more physical as 
it could mimic a disturbance in the gravitational field of 
the star due to the presence of another object. The Zer- 
illi and Newman-Penrose ^4 [2^ gravitational waveforms 
are extracted and compared at different detectors. The 
waveforms damp out on a short timescale. This is con- 
sistent with the expectation that, similar to boson stars, 
soliton stars have only strongly damped modes because 
the scalar field extends to infinity and this allows energy 
to be radiated away rapidly [2^, [Ij, HI] . We also calculate 
the energy radiated in gravitational waves from the Zer- 
illi function. By the end of our simulations the energy as 
a function of time is seen to fiatten out, suggesting that 
the full gravitational waveform has been extracted. 

In Sec. Ill Al we describe the equations for spheri- 
cally symmetric oscillatons. The eigenvalue problem and 
boundary conditions are then discussed in Sec. IIIBI In 
Sec. Ill CI our results for the excited state evolutions are 
presented. Sec. IIII Al discusses the evolution equations 
for the 3D code. The gauge conditions are presented 
next in Sec. IIII 5] with convergence tests and comparison 
to ID results. Sec. IIII CI details the application of small 
nonradial perturbations to ground state soliton stars and 
discusses the extracted gravitational waveforms for dif- 
ferent perturbations. The results are summarized in the 
conclusion. 
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II. SPHERICALLY SYMMETRIC 
CONFIGURATIONS 

A. Mathematical Background 

The action describing a self-gravitating real scalar field 
in a curved spacetime is given by 



IQttG 



where R is the Ricci scalar, g^^ is the metric of the space- 
time, g is the determinant of the metric, $ is the scalar 
field, V its potential of self-interaction, and units with 
h = c = 1 have been used. Greek indices take values 
between and 3 and the Einstein summation convention 
is used. The variation of this action with respect to the 
scalar field leads to the Klein-Gordon equation, which 
can be written as 
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When the variation of Eq. (1) is made with respect to 
the metric g*"", the Einstein's equations G^^ — SttGT^u 
arise, and the resulting stress energy tensor reads 



T^u = - + (3) 

In the present manuscript we focus on the free-field case, 
for which the potential is V{4)) = to^$^, where m is 
interpreted as the mass of the field. 

The spherically symmetric line element is 



ds^ = -N^df + g^dr"^ + r^{de^ 



(4) 



where N{r, t) is the lapse function and g^{r, t) is the ra- 
dial metric function. We use a polar slicing condition and 
this ensures that the line element retains the above form 
throughout the evolution [1, 0]. We follow Refs. 0, [11 
and write the coupled Einstein-Klein-Gordon equations 
as 



AnGrAUJ^^ + + Am^^^) + -(1 - A), (5) 
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where A{r,t) = G(r,t) = [g{r,t)/N{t,r)f and $(r,t) 
is the real scalar field. These equations have no equilib- 
rium solutions with static metric components @. The 
simplest solutions are periodic expansions of the form 



$(t,r) = ^ 02j_i cos((2j ~ l^i), 

jmax 

A{t,r) = ^ A2j(r)cos(2jw0, 
G(i,r) - 



^ G2j(r) cos(2jwi). 



(9) 
(10) 

(11) 



where uj is the fundamental frequency and j'max repre- 
sents the value of j at which the series is truncated for 
numerical computations. The actual solution is an infi- 
nite Fourier expansion of the above form that is conver- 
gent 

For numerical convenience in the code we use dimen- 
sionless variables 



$ ^ \/87rG$, 



r/rn, G Crr? juP' ^ t — s- ujt. 
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B. Boundary Conditions and Eigenvalue Problem 

The system is asymptotically flat: A{r — oo,t) ~ 1 
and <I>(r = oo,t) = 0. Thus the coefficients Ao(oo) = 1, 
Aj (oo) — for j 7^ 0, (pj (oo) = for all j, and Cj (oo) = 
for j ^ 0. The series is truncated at jmax = 2 and the 
coefficients of sin(jt) and cos{jt) are matched in Eqs. ^ 
El). Eq. ((HI) is used as an algebraic relation to determine 
Aq only. The system of equations is solved as an eigen- 
value problem for Cj after specifying the central field 
(/)i(0). The lapse at the boundary gives the value of uj. 
See Ref. [H for further details. One sets t = in Eq. © 
to determine the initial metric components grr{t = 0) = 
A0+A2+A4, gtt{t = 0) = -{Ao+A2+A4)/{Co+C2+C4). 

The configurations we consider are excited states. The 
fields 4>2j-i with j = 1, 2, ... have nodes. The first excited 
state field configurations have 1 node, the second have 2 
and so on. 

After determining the initial configuration, the evolu- 
tion of the system is studied using the ID boson star evo- 
lution code of Refs. [l|,[i^ with the field and its derivative 
having one component (real field) instead of two (com- 
plex field). The mass of the star is determined by the 
value of grr = A{t, r) at the edge of the grid 
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C. Stability of Excited States 

The mass profile of ground state soliton stars have been 
described in Refs. 0,[il. In Fig. [2 we show the mass of 
first excited state oscillatons as a function oi(f)i{0). These 
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FIG. 1: The mass is shown as a function of central density for 
stars in their first excited state. The mass increases mono- 
tonically until it reaches a critical point at (^ic ~ 0.45, and 
Mc « 1.33A'/pi/m after which it starts decreasing. This point 
marks the end of the ^-branch, and the beginning of the f/- 
branch. These stars are more massive than ground state stars. 



stars have two branches similar to the ground state con- 
figurations with the mass increasing to a maximum mass 
of Mc ~ 1.3Afpj/m. This is higher than the maximum 
mass of a ground state star. In general, excited config- 
urations of a given central field density 0i(O) are larger 
and more massive than ground state configurations with 
the same central density. The branch to the left of the 
maximum is traditionally called the S'-branch and the 
branch to the right is called the [/-branch. 

In the case of ground state stars the S-branch is found 
to be stable to radial perturbations Stability here 
means that stars on this branch move to new lower mass 
configurations on the same branch under small pertur- 
bations, [/-branch stars are inherently unstable to small 
perturbations and collapse to black holes. Under large 
perturbations that reduce their mass sufficiently they can 
migrate to the S'-branch Q. 

In this work we study the stability of excited state S'- 
branch oscillaton stars. In Fig. [2] the evolution of the 
radial metric g^r of an S'-branch first excited state star 
is shown. The central density of this configuration is 
(/)i(0) = 0.2828 and no explicit perturbation has been ap- 
plied. The numerical grid (Ar = 0.1) is the only source 
of perturbation. As can be seen from the sharp rise in the 
radial metric (Fig. [SJa)) and collapse of the lapse (Fig. 
[2Kb)) the star is going to form a black hole. The polar 
slicing condition causes the lapse to collapse and the ra- 
dial metric to rise sharply when an apparent horizon is 
approached [l^. Clearly, the star is unstable although 
this configuration lies on the S-branch. 

We next study the time of collapse of one node S- 
branch oscillaton configurations to black holes as a func- 
tion of central density. Excited states are inherently un- 
stable and the discretization error due to the numerical 
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FIG. 2: (a) The lapse and (b) the radial metric are shown 
at different times as a star with a central density of <^i(0) = 
0.2828 collapses to a black hole. The polar slicing condition 
causes the lapse to collapse and the metric to rise sharply as 
an apparent horizon is approached. 



grid is sufficient to make them collapse to black holes or 
migrate to the ground state. To compare the effects of 
radial perturbations on stars with different central den- 
sities, we use the discretization of the grid as a perturba- 
tion with the stars covered by the same numbers of grid 
points. The stars need to be resolved differently because 
of their different central densities. In the following sim- 
ulations the 95% mass radius Ro,^ was considered to be 
the radius of the star. 

The first configuration considered is close to the critical 
mass and has 0i(O) = 0.41 and a radius of i?95 = 14.9/r7i. 
This was covered by 187 grid points with a resolution 
Ar « 0.08. The second configuration has ^i(O) = 0.2828 
star with i?g5 — jm and was covered by the same 
number of grid points with Ar = 0.1. This star col- 
lapsed in a longer time. We took the time of collapse 
(polar slicing does not allow us to find the actual time 
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FIG. 4: (a) The mass of the same (^i(O) = 0.041 star is 
shown as the star migrates to the ground state. The mass 
flattens out as the star is settling towards a final ground 
state configuration. The oscillations in the metric also set- 
tle down to an approximately constant oscillation. The mass 
of the star at the end of the run is M « 0.44Afpi/m, and is 
within 2% of the ground state configuration (0i(O) = 0.11 and 
M — 0.43Mpi/m) that the star is estimated to migrate to. (b) 
The maximum of Qrr is shown as a function of time. The star 
loses more than 30% of its mass during the migration. 



FIG. 3: In (a) the scalar field $ and (b) metric component grr 
are plotted against radius at t = and t « 690, 000/m. The 
initial grr has two maxima (first excited state star; 0i(O) — 
0.041 and Rgs = 51.35/m). By t « 690, 000/m the star has 
migrated to a denser ground state configuration of Rg^ ~ 
16/m with lower mass, (c) The final metric grr is compared 
to that of a ground state t = 0, (^i(O) = 0.11 star. Although 
the ground state grr in slightly lower because the star has not 
fully settled, the overall agreement is very good. 



of formation of a black hole horizon p^]) to be the be- 
ginning of the sharp rise in the metric. Coordinate grid 
stretching [29| (due to the differential rates that test par- 
ticles on geodesic curves fall into the black hole) causes 
the metric function g„ to rise and form a sharp peak. 
The formation of the peak signals the approach of hori- 
zon formation on a short time scale. 

In Table |T] we show the time of collapse to a black hole 
for different S'-branch configurations. The table shows 
the central field density, mass, and radius of each config- 
uration. The time of collapse to a black hole increases 
as the central density decreases until a central density of 
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TABLE I: Physical characteristics of S-branch sohton stars 
in the first excited state are fisted together with the coUapse 
time to black holes. Each star is covered by the same number 
of grid points 187) to facilitate the comparison. The radius 
of the star is taken to be the radius within which 95% of the 
mass of the star is contained. Configurations with a mass 
below 0.84Afp[/m lose enough mass through scalar radiation 
to migrate to the ground state. 
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FIG. 5: (a) The radial metric Qrr of a migrating five node 
soliton star configuration (0i(O) = 0.006; M = 0.79MpJm) 
is displayed as a function of radius on two timeslices. The 
t = radial metric has 6 maxima, which is characteristic to 
5th excited state configurations. By t ~ 172,000/m the star 
has migrated to the ground state. Although, the migration 
process is not yet complete by this stage, the star has lost 29% 
of its mass through scalar radiation and the radial metric has 
only one maximum, (b) The central density of the same star is 
shown at t — and in a intermediate state at t ~ 102, 000/m. 
In this intermediate state the star has four nodes (central 
density has five maxima). We plot the density of a four node 
star with <^i(0) — 0.007 for comparison. 



about 01 (0) — 0.075. For central densities of 0.075 and 
lower the stars did not collapse but migrated to config- 
urations in the ground state. Although a first excited 
state star of central density = 0.075 is more massive 
than any ground state star (maximum mass of ground 
state configurations is about 0.61Mp^/m) it lost the ex- 
cess mass through scalar radiation and settled in the 
ground state. Excited state S'-branch stars are clearly 
inherently unstable and either collapse to black holes or 
migrate to the ground state. 

In Fig. [3] and Fig. 2] we show a first excited state star 



with 4>i{0) ~ 0.041 of mass 0.655Af|[/r7i migrating to 
the ground state. No explicit perturbation was applied. 
In Fig. El^a) the initial (one node) and final (no nodes) 
field configurations are displayed. In Fig. [3{b) the initial 
and final radial metric shown. The initial radial 

metric has two maxima, which is characteristic of a first 
excited state configuration. The final radial metric {t = 
218,9987r/m) has one maximum indicating the star has 
gone to the ground state. The final configuration shown 
corresponds to a time which is an integral multiple of tt. 
This is to ensure that we get the same phase as the t = 
configuration to facilitate comparison. In Fig. [3jc) the 
radial metric of the star at the end of the run is plotted 
together with the radial metric of a ground state star 
of central density 0i(O) = 0.11. The final radial metric 
of the migrating star is very close, to this ground state 
configuration. For comparison the mass of the star at the 
end of the run was about M — OAAAIp^/m as compared 
to a mass of M = 0.43M|i/to for a 0i(O) = 0.11 star. 
The migrating star is settling to a configuration close to 
the (j)i{0) = 0.11 ground state star. 

The mass loss is calculated using Eq. (IT^ for the above 
configuration and is plotted as a function of time in Fig. 
|4fa). The star loses more than 30% of its mass through 
scalar radiation. The amount of mass loss decreases in 
time as the star settles to its final ground state configu- 
ration. In Fig. SKb) the maximum of the metric function 
Qrr of a perturbed excited state star is plotted against 
time. As the star settles into its new ground state config- 
uration, it oscillates at two different frequencies. One is 
approximately the fundamantal unperturbed oscillation 
of period tt/uj (icodc = tt). The second has much higher 
period that changes as the star goes through different in- 
termediate configurations and eventually settles into the 
quasinormal mode of the ground state configuration Q . 

We next simulate the migration of an oscillaton in 
the 5th excited state to the ground state. The ini- 
tial configuration is that of a five node excited state S'- 
branch star of central density (/)i(0) = 0.006 and mass 
M = 0.792M|[/m. As usual, no explicit has been applied 
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perturbation other than that due to the discretization of 
the numerical grid (resolution Ar = 0.1). In Fig. [5l^a) 
we show the radial metric grr of the star at t = 0, and 
t « 172, 000/m. The initial configuration (< = 0) has six 
maxima in the radial metric. By the end of the simulation 
the radial metric has just one maximum indicating the 
star has gone to the ground state (the mass at this stage 
is M = 0.463Mp[/m.) During the evolution the star ap- 
pears to pass through an intermediate four node state. 
This can be seen in Fig.[5]^b). In the figure p/ cos{2t) [p 
is the density of the star and the division by cos(2i) is in 
order to facilitate comparison as it helps avoid phase dif- 
ferences) is plotted aX t — and t k, 102, 000/m against 
the t = Q density of a four node (t>i{0) — 0.007 star. The 
evolving star initially had a density p with six maxima 
characteristic of a five node star. By time t = 102, 000/m 
it has lost one of these maxima and appears to be in an in- 
termediate four node state. We have plotted the density 
of a four node (t>i{Q) — 0.007 star to compare the posi- 
tions of the peaks and profile of p. The positions of the 
peaks are close to each other (especially the first three) 
but the sizes of the peaks are not. The star probably 
cascades through superpositions of lower excited config- 
urations due to the inherent non-linearity of the system. 



III. NUMERICAL SIMULATIONS OF GROUND 
STATE SOLITON STARS IN 3D 

A. 3D Evolution equations 

In order to find solutions to the Einstein-Klcin-Gordon 
system of equations we use the 34-1 decomposition of 
Einstein's equations, for which the line element can be 
written as 



a'dt^ + -labidx"- + I3''dt){dx'' + fi'dt), (14) 



where jab is the three dimensional metric; from now on 
latin indices label the three spatial coordinates. The 
functions a and (3°" in Eq. (fT4|) are freely specifiable 
gauge parameters, known as the lapse function and the 
shift vector respectively, and 7 is the determinant of the 
3-mentric. Throughout the rest of the paper the standard 
general relativity notation is used. The Greek indices run 
from to 3 and the Latin indices run from 1 to 3. In this 
section we use geometric units with G = c = 1. 

The Klein-Gordon equation can be written as a first- 
order evolution system by defining two new variables 
in terms of combinations of their derivatives: tt — 
{y/j/a){dt(f> — P'^dccj)) and ipa = da(j>- With this nota- 
tion the evolution equations become 



J2 



(15) 



1 dV 

On the other hand, the geometry of the spacetime is 
evolved using the BSSN formulation of the 3-1-1 decom- 
position. According to this formulation, the variables to 
be evolved are 5* 



-4* 



{Kab 



. In 7/12, 7a6 = e-4*7,b, K = ^'"'Kab, 
labK/3) and the contracted Christof- 



fel symbols = 7 F^^, instead of the usual ADM vari- 
ables Jab and Kab- The evolutio n eq uations for these new 
variables are described in Refs. 123. [2311: 
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dtK 
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(18) 
(19) 
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(/3=r'c-2r''r?c + 3r'/3%),(20) 



d_ 

"a? V" ' - ' ^= ■ 3 

where Da is the covariant derivative in the spatial hyper- 
surface, T is the trace of the stress-energy tensor ^ and 
the label TF indicates the trace-free part of the quantity 
in brackets. The coupling between the evolution of the 
BSSN variables and the variables describing the evolu- 
tion of the scalar field is first order. That is, Eqs. (5) 
are solved using the method of lines with a modified ver- 
sion of the second order iterative Crank-Nicholson (ICN) 
integrator (see Ref. H^l). After a full time step the stress- 
energy tensor in Eq. (3) is calculated and used to solve 
the BSSN evolution equations with an independent evo- 
lution loop based on the standard second order ICN [3l[ . 

In order to set up the correct scaled quantities to 
be evolved we use dimensionless variables. For an 
equilibrium configuration the real scalar field can be 
expanded in a Fourier series of the form $(r, t) — 
X]j=T "^2^-1 (»') cos((2j — 1)^^), where the total central 

density is (f){Q) = X]j=r '/'2j-i(0). This implies that the 
stress energy components given in Eq. (3) also have a 
periodic time dependence. The characteristic frequency 
Lo together with the mass of the boson m provide the 
natural dimensionless units, 



''code = mr/Mpi, tcodo = ujt, 



(21) 



47r 



^codc 



Mpi ^' 



^codc ^ 



which are the ones used within the present analysis. For 
further code details refer to Ref. (2l| . 



B. Slicing Condition and Convergence Tests 

In this section we introduce the shcing condition used 
for our 3D simulations. We base it on the observed dy- 
namics of the trace of the extrinsic curvature tensor K for 
oscillaton eigenstates, and find that the shcing is accu- 
rate for evolving spherically symmetric oscillatons with 
small nonradial perturbations. For soliton star eigen- 
states K{t = 0) = Krr{t = 0) = 0. However, as time 
evolves the metric function g^r changes with time and 
consequently K oscillates periodically, given by 



K = 



-2q -2a 
We denote K expanded to order j 



= Y,K2jSm{2jt). (22) 



as 



^0„,ax) ^ K2j sin{2 jt). 



(23) 



In the case of spherical symmetry with the gauge condi- 
tion provided by polar slicing one can determine K2j by 
using the Fourier expansion of and a = ^ AjC from 
Eqs. (jlOjllip in the expression for K (Eq. (HH)). One can 
match the coefficients of sin(2jf) for each j and obtain 
in terms of a combination of A2^ and C2, . The first 
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two coefficients K2 and K/^ are approximately given by 
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(24) 



(25) 



For 3D simulations we propose to enforce that K take 
the value ^(J— ) of Eg. and do so in stable man- 

ner. Various authors [32j have noted that one can enforce 
maximal slicing K = = K in sl stable manner by solving 



dK 
'dt 



= -cK, 



(26) 



where c is a positive constant. Analogously, we drive K 
to the expanded iC^-'"'"'^ for an eigenstate using 

l(^-*-"" 

Using the identity 



M = -c K - K'^^" 



(27) 



V^a - KabK^^a - ^it{S + p)a 



dK 



(28) 



where p is the energy density and S is the trace of the 
spatial stress tensor, we rewrite Eq. ((27)) as an elliptic 
equation to solve for the lapse 



V^a - i^a&i^ a - 47r(S' + p)a = - 



dt 

c(k - K^^" 



(29) 




FIG. 6: The extrinsic curvature K is shown for a Model 5*1 
star on different timeslices as a function of r. The solid lines 
represent the results obtained from a ID spherically sym- 
metric simulation, while the symbols (filled and empty ci- 
cles, squares and triangles) are 3D data on the same times- 
lices from a Model Si star simulation on a 192^ grid with 
Aa; = Ay — Az = 0.15 resolution. Good agreement is ob- 
served. 



We refer to the slicing obtained from this gauge condition 
(Eqs. (P7)) and (PP]) ) as truncated Fourier slicing. 

For the 3D simulations described in this section and 
the subsequent section, we use a ground state oscillating 
soliton star with 



61(0) = 0.20, Lu = 0.91 



M 



0.57^, 



R 



'95 



Mil 



(30) 



and we label this as the Model Si oscillaton for future 
reference. For initial data within the 3D simulations we 
take high resolution ID spherical data for this model and 
interpolate it onto the 3D Cartesian grid using bspline 
interpolation. All our 3D simulations are performed in 
octant symmetry. 

Fig. [6] shows if as a function of r on different time 
slices for a Model 5*1 soliton star simulated on a 192'^ 
three-dimensional grid with a resolution of Ax — Ay — 
Az ~ 0.15. This is compared with the ID result for 
a Ar = 0.01 resolution. The spherically symmetric ID 
evolution of K for a soliton star eigenstate is a periodic 
oscillation between an upper bounding profile at approxi- 
mately t « 7r/4 and a lower bounding profile at t w 37r/4. 
The oscillations of K are observed to match the expected 
time dependence of Eq. ([^5]) and the agreement between 
the ID and 3D results is good, indicating the truncated 
Fourier slicing condition is enforced accurately. 

In order to demonstrate the advantages of applying 
this gauge condition, we simulate the same Model Si 
star on a 96'^ grid with Ax = Ay = Az = 0.30 with 
truncated Fourier slicing and two widely used slicing 
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FIG. 7: The L2-norm of the Hamiltonian is displayed as a 
function of time for three slicing conditions: (1) 1+Iog slicing, 
(2) Maximal Slicing and (3) Truncated Fourier Slicing. The 
latter is clearly smaller and leads to more stable simulations. 



conditions: maximal slicing (K = 0) and 1+log slicing 
{dt(x — —2aK). Fig. [7] compares the L2-norm of the 
hamiltonian constraint for these three simulations. It can 
be seen that the Hamiltonian for the truncated Fourier 
slicing is by far the lowest (smaller by approximately a 
factor of 2 than 1+log and Maximal slicing) and behaves 
is a stable manner. These three test runs are carried 
out to a time of 47r; for longer simulations that study 
perturbed soliton stars we find that the use of the trun- 
cated Fourier slicing condition is necessary to control the 
coordinate drifting of metric functions to an acceptable 
tolerance and to perform the gravitational wave extrac- 
tion accurately. 

A further test of the accuracy of the evolution code 
using the truncated Fourier slicing condition is the con- 
vergence test shown in Fig. [51 Fig. [SJa) displays the 
first oscillation of the maximum of the density p of the 
Model 5*1 star converging to the ID result as the 3D 
resolution is improved. We take three different resolu- 
tions: 48^ grid with Aa; = Ay = Az = 0.60 resolution, 
96^ grid with Ax = Ay = Az = 0.30 and 192^ grid 
with Ax = Ay = Az = 0.15 resolution. These are com- 
pared to a high resolution spherically symmetric simu- 
lation (Ar = 0.01) with the ID code. Fig. [5]^b) shows 
the convergence of the hamiltonian constraint, which is 
observed to be approximately second order. 



C. Nonradial Perturbations of Soliton Stars and 
Waveform Extraction 

The focus of this section is the study of small non- 
spherical perturbations applied to spherically symmetric 
soliton stars and the extraction of gravitational wave- 
forms for such systems. We first perturb the star by 
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FIG. 8: (a) The maximum density p of the Model Si soli- 
ton star is shown converging to the spherically symmetric ID 
result, (b) The L2-norm of the Hamiltonian constraint is dis- 
played for the three different resolutions. The code is observed 
to exhibit roughly second order convergence. 



imposing an asymmetry in the grid resolution. We simu- 
late the Model Si star using a different resolution in the 
z-direction than in the x and y directions, while keeping 
the distance from the origin to the boundary on each axis 
the same. Fig. [5] shows the Zerilli waveform extracted at 
r = 43Mpj/m for three different resolutions labeled as 
(1) high resolution for Ax = Ay = 0.196, Az = 0.201 
with = Hy = 244, Uz = 238, (2) medium res- 
olution for Ax = Ay = 0.294, Az = 0.3015 with 
TT-x = ny = 164, Uz = 160, and (3) low resolution 
for Ax ^ Ay = 0.392, Az = 0.402 with Ux ^ Uy ^ 
122, Uz — 119. The waveforms are observed to have the 
same frequency but different amplitudes. For this series 
of simulations the uneven resolution is itself the perturba- 
tion, and thus the waveform should converge away as the 
resolution is improved. We observe the appropriate sec- 
ond order converge rate, with a factor of four difference 
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FIG. 9: The i = 2,m = Zerilli function extracted at r = 43 
is shown for a perturbed Model Si sohtion star at different 
resolutions labeled as: High for Aa; — Ay = 0.196, Az — 
0.201 with = riy = 244, = 238, Medium for Ax = 
Ay = 0.294, Az = 0.3015 with = = 164, = 160, 
and low resolution for Ax = Ay = 0.392, Az = 0.402 with 
tlx = riy = 122, — 119. The perturbation in this case is the 
asymmetric resolution itself and amplitude of the waveform 
is seen to converge away. The convergence rate is roughly 
second order. 



between the waveform amplitudes of the high resolution 
and the low resolution case. Each of the simulations is 
performed using the truncated Fourier slicing condition 
and is run up to a time t « 60 . 

We next study an oscillating soliton star with a per- 
turbation applied to the radial metric component g^r of 
the form 

Sgrrir, 9, <P) - e2^f(r/R)Y2o{e, ^)grr{r) (31) 

where £20 is a constant much less than unity and the 
function f{r/R) = (r/Rp)'^ exp{—cr'^). The constants c 
and Rp are chosen to localize the perturbation on the 
star and also not perturb the metric functions at the ori- 
gin. This perturbation does not significantly alter the 
mass of the star. We expect that physical perturbations 
of soliton stars can be easily described as linear combina- 
tions of spherical harmonics. This type of perturbation 
in the metric may mimic the effect that gravitation of 
another star has on the oscillaton. In the current study 
we simply apply the perturbation of Eq. (PT|) without re- 
solving the initial data problem. This procedure causes 
the constraints to no longer be obeyed and hence makes 
the perturbation nonphysical. A more complete treat- 
ment that incorporates an initial value problem solver 
that couples a scalar field and metric perturbation is the 
subject of future work. 

We study the effects of a perturbation with £20 = 10"**, 
c = 1/16 and Rp = 5 on a Model 5*1 soliton star. In Fig. 
[TUf a) the £ = 2,m = Zerilli functions extracted at 
a detector radius of r = A^Mp^/m are shown for three 




t 

FIG. 10; (a) The ^ = 2, m = Zerilh function extracted at 
r = 45 is shown for a Model Si solition star under a Y2Q 
perturbation at three different resolutions: High for Aa:: = 
Ay — Az — 0.196 with n^c = riy = riz = 244, Medium for 
Aa; — Ay = Az = 0.294 with ~ riy = — 164, and low 
resolution for Aa; = Ay = Az = 0.392 with rij, = = = 
122. The waveforms have roughly the same frequency and 
amplitude, (b) The Zerilli energy is displayed and is observed 
to flatted out on the same timescale as the gravitational wave 
signal. 



different resolutions. The waveforms oscillate through 
a couple of nodes and then damp to zero on a fairly 
short timescale. This is similar to the case of boson stars 
[ifi, 26] and consistent with the expectation that a star 
with a scalar field exten ding to infinity would have only 
rapidly damped modes [la [13, The energy emit- 

ted in gravitational waves (see Fig. [TUTb)) is calculated 
and found to fiatten out at about 5 x 10~^Mp[/TO. This 
suggests that the whole gravitational waveform has been 
extracted. 

We now investigate gravitational waveforms using the 
Newman-Penrose scalar ^1^4 for the simulation described 
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FIG. 11: The Newman-Penrose scalar '^4 multiplied by the 
radius at the detector is shown at three different detectors. 
The signals have been timeshifted by the distance in flat space 
between detectors to match the last detector. The agreement 
in both frequency and amplitude is good. The inset shows 
that the high frequency oscillations are noise and converge 
away with resolution. 



above. Fig. [TT] shows r^4 at three different detectors. 
The waveforms are timeshifted by the distance in flat 
space between each detector and the last detector to fa- 
cilitate comparison. The frequencies and amphtudes of 
the waveforms at the three detectors are in good agree- 
ment. This is consistent with the expected behavior of 
^4 oc 1/r in the far wave zone (25| . High frequency 
oscillations appear at late times. They are likely an arti- 
fact of the numerical discretization. The frequency and 
amplitude of these oscillations arc resolution dependent 
(see inset); these data show the amplitude to converge 
away to second order. A more advanced technique such 
as adaptive mesh refinement that provides higher resolu- 
tion would enable us to extend the convergence test to 
finer grid resolutions. 

Within all our 3D simulations we use the L2-norm of 
the Hamiltonian constrant as a measure of the error. This 
is typically of the order of a few x 10~^ and below 1 x 10^^ 
for all the simulations in this section. 



IV. CONCLUSIONS 

For the first time excited state soliton stars have been 
studied. Excited state configurations have field configu- 
rations characterized by nodes. For a given central field 
density (f>i (0) , these stars are typically larger (greater ra- 
dius) and more massive than corresponding ground state 
stars. We find that their mass profile is similar to that 
of ground state configurations; they have two branches 
separated at the maximum mass. The branch to the 
left of this maximum is traditionally called the ^-branch 
because ground state configurations on this branch are 



stable to radial perturbations in the sense that they mi- 
grate to new configurations on the same branch. In spite 
of this similarity in the mass profile, we find that all 
excited state configurations arc inherently unstable. Un- 
like ground state configurations, excite state stars on the 
5-branch do not migrate to new configurations on the 
excited state 5-branch when perturbed. They can either 
migrate to the ground state if they lose enough mass 
through scalar radiation or collapse to black holes. 

Higher excited state configurations migrate to the 
ground state sometimes cascading through intermediate 
excited states. We show one such migration for a five- 
node excited state star and find the four-node intermedi- 
ate state that it is nearest to during the migration pro- 
cess. Studying the stability of excited states is very im- 
portant because they may be intermediate states during 
the formation process of soliton stars. 

Also, for the first time, the numerical evolution of 
ground state oscillaton stars is conducted on a full 3D 
grid, allowing the study of nonradial perturbations with 
gravitation waveform extraction. These simulations are 
challenging because soliton stars are very dynamic, with 
no equilibrium configurations having static metric com- 
ponents. In this paper we compare slicing conditions pre- 
viously used for 3D boson star simulations with a new 
slicing condition specifically designed for soliton stars. 
We find that the latter is more accurate and more stable. 
We see the energy density in 3D converging with reso- 
lution to that in ID for the spherically symmetric case 
and the highest resolution 3D simulation matching the 
ID result very accurately. 

In 3D we explore two types of nonradial axisymmetric 
perturbations. First, the spherically symmetric ID eigen- 
state data is interpolated on the 3D grid with an imposed 
asymmetry in the grid resolution Ax — Ay 7^ Az. No 
explicit perturbation is applied. This has the advantage 
that the initial data automatically satisfies the Einstein- 
Klein-Gordon equations. However, this perturbation is 
not easy to interpret physically. We successfully extract 
the Zerilli gravitational waveforms for a series of simula- 
tions at different resolutions and observe the amplitude 
of the signal due to this resolution-based perturbation 
to converge away with improved resolution. The con- 
vergence rate is approximately second order as expected. 
A second perturbation is applied explicitly to the met- 
ric function g„ and is proportional to the Y2Q spheri- 
cal harmonic. Both the ^ = 2, m = Zerilli function 
and the Newman Penrose scalar ^['4 are extracted. This 
perturbation more closely represents a physical distur- 
bance of the star and leads to a finite signal that does 
not converge away as in the previous case. The wave- 
forms damp rapidly and thus the complete gravitational 
waveform could be extracted on a relatively short time 
scale. This is consistent with expectations from pertur- 
bation theory for stars that have a scalar field extending 
to infinity such as soliton stars and boson stars [1^. We 
also measure the energy radiated in gravitational waves 
as a function of time and find that it fiattens as the signal 
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damps out. 

This is an exploratory investigation that is meant as 
a starting point for future work. Possible improvements 
include using an adaptive mesh refinement technique for 
better resolution and finding a more general Gauge con- 
dition that would be applicable to strongly perturbed 
soliton stars. These would enable the study of more ad- 
vanced scenarios involving real fields such as soliton star 
collisions. 
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